StatQuest生物统计学专题 - PCA的奇异值分解过程

特征值分解（或奇异值分解）的原理PCA是如何同奇异值分解联系在一起的PCA的奇异值分解的R过程

特征值分解（或奇异值分解）的原理

特征值分解和奇异值分解两者有着很紧密的关系，特征值分解和奇异值分解的目的都是一样，就是提取出一个矩阵最重要的特征。

特征值分解

如果说一个向量v是方阵A的特征向量，将一定可以表示成下面的形式这时候λ就被称为特征向量v对应的特征值，一个矩阵的一组特征向量是一组正交向量。特征值分解是将一个矩阵分解成下面的形式： 

其中Q是这个矩阵A的特征向量组成的矩阵，Σ是一个对角阵，每一个对角线上的元素就是一个特征值。

奇异值分解

下面重点谈谈奇异值分解。特征值分解是一个提取矩阵特征很不错的方法，但是它只是对方阵而言的。

奇异值分解是一个能适用于任意的矩阵的一种分解的方法：

假设A是一个M * N的矩阵，那么得到的U是一个M * M的方阵（里面的向量是正交的，U里面的向量称为左奇异向量），Σ是一个M * N的矩阵（除了对角线的元素都是0，对角线上的元素称为奇异值），V’(V的转置)是一个N * N的矩阵，里面的向量也是正交的，V里面的向量称为右奇异向量） 

那么奇异值和特征值是怎么对应起来的呢？首先，我们将一个矩阵A的转置 左乘A，将会得到一个方阵，我们用这个方阵求特征值可以得到： 

这里得到的v，就是我们上面的右奇异向量。此外我们还可以得到： 

这里的σ就是上面说的奇异值，u就是上面说的左奇异向量。

来源于引文2：奇异值分解(SVD)详解及其应用。

PCA是如何同奇异值分解联系在一起的

在上一节的分享StatQuest生物统计学专题 - PCA中，我们知道PCA的主成分求解的直观理解就是依次对数据进行拟合：

PC1是拟合一条过原点的最佳直线；

PC2是在同PC1正交的方向上拟合另一条最佳直线；

PC3是在同PC1、PC2平面垂直的方向上拟合一条最佳直线；

PC4是在同PC1、PC2、PC3相正交的方向上拟合。

以此类推，最多有n个主成分（n为样本数量）。

如下图的求解PC1为例，先看单个点的情况，那么最佳拟合直线（红色虚线）的指标就是使得点与直线之间的距离a最小，由于斜边c固定不变，所以也就是使得距离b最大。

所以考虑全部点之后，最佳拟合直线就是使得全部点的b之和最大，以b平方表示距离，也就是使得最大。

在PCA中，最佳拟合直线的就是特征值，就是奇异值。

而最佳拟合直线所代表的向量就是特征向量或奇异值向量。

如上一节的例子中："PC1由4份Gene1和1份Gene2构成" 。向[4,1]其实类似于奇异值向量，只不过奇异值分解时，向量是标准化的，所以准确来说向量[0.97,0.242]才是PC1的奇异值向量。

PCA的奇异值分解的R过程

在上一节中，已经知道R中是使用prcomp函数进行PCA分析。

#构造数据
data.matrix <- matrix(nrow=100, ncol=10)
colnames(data.matrix) <- c(paste0("wt", 1:5),  paste0("ko", 1:5))
rownames(data.matrix) <- paste0("gene", 1:100)
for (i in 1:100) {
  wt.values <- rpois(5, lambda=sample(x=10:1000, size=1))
  ko.values <- rpois(5, lambda=sample(x=10:1000, size=1))

  data.matrix[i,] <- c(wt.values, ko.values)
}
data.matrix <- t(data.matrix)
#pca分析
pca <- prcomp(data.matrix, scale=TRUE)
#pca中有三个参数x、sdev、rotataion
pca$x #各个主成分
pca$sdev #主成分的标准差
pca$rotation #主成分的载入因子

下面将使用奇异值分解来完成PCA，奇异值分解函数使用svd来实现：

svd.stuff <- svd(scale(t(data.matrix), center=TRUE))

主成分

prcomp函数返回对象pca，pca中的x就是主成分。

svd函数返回对象svd.stuff，svd.stuff的u就是主成分，由于奇异值分解进行了标准化，因此再乘以svd.stuff的d就是相应的主成分。

如下图，两者的PCA图完全相同：

opar <- par(no.readonly=T)
par(mfrow=c(1,2))
plot(pca$x[,1],pca$x[,2],xlab='PC1',ylab='PC2',main='Principal Component Ananlysis')
plot(svd.stuff$u[,1]*svd.stuff$d[1],svd.stuff$u[,2]*svd.stuff$d[2],xlab='PC1',ylab='PC2',main='Singular Value Decomposition')
par(opar)

载入因子

载入因子就是各个主成分的成分构成，也就是奇异向量。

pca对象中rotation代表了载入因子，共有10个主成分，每个主成分是100个gene值的线性组合，这些线性组合的系数就是载入因子，就是奇异向量。

svd.stuff中v就是载入因子。

如下图，pca的PC1载入因子同svd的PC1载入因子完全相同。

plot(pca$totation[,1],svd.stuff$v[,1],xlab='pca',ylab='svd')

方差

svd.stuff中有一个参数d，d就是每一个主成分所能解释的标准差，也就是奇异值（d的平方，也就是特征值）。由于每一个主成分就是最佳拟合直线，所以这个标准差其实就是上面提到的投影点到原点距离平方和的平方根。

pca中有一个参数为sdev，这个是每个主成分的标准差，和svd.stuff$d的意义是一样的。不同的地方在于，pca对sdev进行了无偏估计，而svd.stuff的d是直接进行开方所得：

所以，将svd.stuff的d进行n-1修正后，就会和pca$sdev完全相同。

svd.stuff.sdev <- sqrt(svd.stuff$d^2/(ncol(data.matrix)-1))
plot(pca$sdev,svd.stuff.sdev,xlab='pca',ylab='svd')

专题以往文章

参考资料

StatQuest课程：https://statquest.org/video-index/

奇异值分解(SVD)详解及其应用：https://blog.csdn.net/shenziheng1/article/details/52916278

---

 猜你喜欢

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。

▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。